#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
# Simulation functions 
#   SIMULATION functions
#          1)      mrw.sim  : mu random walk 1
#          2)     mrwd.sim  : mu random walk 1 with a drift
#          3)    mrwrd.sim  : mu random walk with a random drift 1
#          4)   mrwAll.sim  : mu random walks with different orders
#          5)      mar.sim  : mu AR process (stationary)
#          6)      mAR.sim  : mu AR process (non stationary process are allowed)
#          7)      yrw.sim  : random walk for y
#          8)     yrwd.sim  : random walk with a drift for y
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------
#----------------------------------------------------------------------------------------  
# functions to generate y values from random walk senarios
#-------------------------------------------------------------------
# simularing random processes 
#---------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------
#---------------------------------------------------------------------------------------
# function 1
# RANDOM WALK for mu order 1
# note that the random walk is on mu not y
#-------------------------------------------------------------------
mrw.sim  <- function(N=100, mu=0, sige=10, sigb=1, plot=FALSE)
{
  wn1 <- ts(rnorm(N, 0, sigb))
    T <- ts(mu+cumsum(wn1))
    y <- ts(T+rnorm(N, m=0,s=sige))
   if (plot) {plot(y); lines(T, col="red")}               
    attr(y, "mean") <- T
    y
}
#-------------------------------------------------------------------
 mrw.plot <- function(N=100, k=5, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- mrw.sim(N, ...)
matplot(y, type="l")   
}
#------------------------------------------------------------------
#------------------------------------------------------------------
#------------------------------------------------------------------
# function 2 
# this is random walk for mu with a drift
mrwd.sim <- function(N=1000, mu=100, d=.1, sige=10, sigb=1, plot=FALSE)
{
  wn1 <- ts(rnorm(N, 0, sigb))
    D <- mu+d*(1:N)
    y <- ts(D+cumsum(wn1))+rnorm(N, m=0,s=sige)
  if (plot)
  {
    plot.ts(y)
  lines(D, col="red")
  }
 attr(y, "drift") <- D
    y  
}
#------------------------------------------------------------------
mrwd.plot <- function(N=100, k=5, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- mrwd.sim(N, ...)
matplot(y, type="l")   
}  
#------------------------------------------------------------------
#------------------------------------------------------------------
#------------------------------------------------------------------
# function 3
# this is random walk for mu with a random drift
mrwrd.sim <- function(N=100, mu=0, d=.5, sige=10, sigb=1, sigd=1, plot=FALSE)
{
  no1 <- ts(rnorm(N, 0, sigb))
  no2 <- ts(rnorm(N, 0, sigd))
    D <- ts(d+cumsum(no2)) 
    #M <- filter(Z, 1, "r", init=mu)
     M <- mu+D+cumsum(no1)
    y <- ts(M+rnorm(N, m=0,s=sige))
  if (plot)
  {
    plot.ts(y)
  lines(T, col="green")    
  lines(I(M), col="red")
  lines(D, col="blue")  
  }
  attr(y, "mean") <- M
    y  
}
#-------------------------------------------------------------------
  mrwrd.plot <- function(N=100, k=5, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- mrwrd.sim(N, ...)
matplot(y, type="l")   
}
#-------------------------------------------------------------------
#-------------------------------------------------------------------
#-------------------------------------------------------------------
# function 4
# RANDOM WALK for mu of different  order
 mrwAll.sim <- function(N, mu0 = 0, sige = 1, sigb = 1, order=1, plot=FALSE)  
 {
# local function to obtain the coef   
    f <- function(k)
    {   
       coef <- rep(0,k+1)
       for (j in 0:k) # this is based on the coefficient  of Diff oparator Rue and Held (2005) p 87
       coef[j+1] <- (-1)^k *(-1)^j * choose(k,j)
      -coef[1]*coef[-1] # this takes the right coef
     } 
#-----
   require(gamlss)
    if (order > 6||order<0 ) stop("order should be between 1 to 6")# you dont realy need this
     N1 <- N+order  # need N+length(phi) x's
 # generate mu(t)= phi*mu(t-1)+e(t), t=1,...,N
 #  and x(0) from the stationary distribution:
    U <- rnorm(N1,0,sigb)
    #switch(order, 
    #       phi <- 1, 
    #       phi <- c(2, -1),
    #       phi <- c(3, -3, 1))
     phi <- f(order)
      mu <- filter(U, phi, "r")
       T <- mu0+mu[-c(1:order)]
    v <- rnorm(N,0,sige)         # obs noise    
    y <- ts(T+v)
    if (plot) { plot(y); lines(T, col="red") }              
    attr(y, "mean") <- T
    y   
 }
#-------------------------------------------------------------------
 mrwAll.plot <- function(N=100, k=5, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- mrwAll.sim(N, ...)
matplot(y, type="l")   
}

#-------------------------------------------------------------------
#-------------------------------------------------------------------
#------------------------------------------------------------------
# is probaly is not correct
#simRWAll <- function(N, init=0, sige=1, sigb=1, order=1, plot=FALSE)
# {
#  require(gamlss)
#     beta. <- rNO(N-order, mu=0, sigma=sigb)
#     if (length(init)!=order) init <- rep(init, order)
#     gamma. <- diffinv(beta., differences=order, xi=init) 
#    #gamma. <- ts(c(init, init+cumsum(beta.)))
#    #cbind(gamma., gamma1)
#    y <- ts(gamma.+rNO(N, mu=0,sigma=sige))
#     if (plot) plot(y)
#    y
#  }
#--------------------------------------------------------------------
#--------------------------------------------------------------------
#-------------------------------------------------------------------
# function 5
 # AR for mu ALL models
# only stationary AR are produce here
 mar.sim <- function(N, mu0=0, sige=1, sigb=1, phi=c(.5), plot=FALSE)  
 {
   require(gamlss)
   lphi <-length(phi)
     N1 <- N+lphi  # need N+length(phi) x'sasssa
 # generate mu(t)= phi*mu(t-1)+e(t), t=1,...,N
 #  and x(0) from the stationary distribution:
   mu <- arima.sim(n=N1, list(ar = phi, sd=sigb), )     # here you have x0 to x100
    v <- rnorm(N,0,sige)         # obs noise    
    y <- ts(mu0+mu[-c(1:lphi)]+v)
    if (plot) plot(y)
    y   
 }
#--------------------------------------------------------------------   
#-------------------------------------------------------------------
 mar.plot <- function(N=1000, k=5, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- mar.sim(N, ...)
matplot(y, type="l")   
}

#--------------------------------------------------------------------
#--------------------------------------------------------------------
#--------------------------------------------------------------------
# function 6
 # AR for mu ALL models
# non stationary AR are produce here
 mAR.sim <- function(N, mu0=0, sige=1, sigb=1, phi=c(.5), plot=FALSE)  
 {
   require(gamlss)
   lphi <-length(phi)
     N1 <- N+lphi  # need N+length(phi) x's
 # generate mu(t)= phi*mu(t-1)+e(t), t=1,...,N
 #  and x(0) from the stationary distribution:
    U <- rnorm(N1,0,sigb)
   mu <- filter(U, phi, "r")
    v <- rnorm(N,0,sige)         # obs noise    
    y <- ts(mu0+mu[-c(1:lphi)]+v)
    if (plot) plot(y)               
    y   
 }
#-------------------------------------------------------------------
 mAR.plot <- function(N=1000, k=5, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- mAR.sim(N, ...)
matplot(y, type="l")   
}
#-------------------------------------------------------------------
#-------------------------------------------------------------------
#------------------------------------------------------------------
# function 7
yrw.sim  <- function(N=1000, mu=0, sigma=1,  plot=FALSE)
{
  wn1 <- ts(rnorm(N, 0, sigma))
    y <- ts(mu+cumsum(wn1))
 mean <- rep(mu, N) 
   if (plot) {plot(y); lines(mean, col="red")}               
    #attr(y, "mean") <- T
    y
} 
#------------------------------------------------------------------
yrw.plot <- function(N=100, k=10, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- yrw.sim(N, ...)
matplot(y, type="l")   
}
#-------------------------------------------------------------------
#-------------------------------------------------------------------
#------------------------------------------------------------------
# function 8
yrwd.sim <- function(N=1000, mu=0, d=.1, sigma=1,  plot=FALSE)
{
  wn1 <- ts(rnorm(N, 0, sigma))
    D <- mu+d*(1:N)
    y <- ts(D+cumsum(wn1))
  if (plot)
  {
    plot.ts(y)
  lines(D, col="red")
  }
 attr(y, "drift") <- D
    y  
}
#----------------------------------------------------
yrwd.plot <- function(N=100, k=5, ...)    
{
  y <- matrix(0, nrow=N, ncol=k)
for (i in 1:k) y[,i] <- yrwd.sim(N, ...)
matplot(y, type="l")   
}
